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PREFACE 

The  research  described  in  this  report  represents  the  con¬ 
tinuation  of  an  extended  investigation  of  numerical  errors  in 
dynamical  veather  prediction,  initiated  and  performed  under  the 
previous  contracts  AF  19(6o4)-3886  and  A F  19(6o4)-4965  during 
the  period  1  July  1958-30  June  i960. 
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Miss  L.  Fujinaka,  13  June  i960  -  3°  September  i960 
Mr.  J.  R.  Hammond,  9  September  i960  -  15  June  1962 
Mr.  E.  G.  Lucero,  25  January  1962  -  25  May  1962 
Miss  M.  H.  Oberlinder,  1  July  i960  -  31  December  1962 
Mr.  C.  A.  Riegel,  1  July  i960  -  30  September  i960 
Mr.  P.  J.  Wlnnlngboff,  9  September  i960  -  30  June  1961 

Bnployed  on  an  occasional  hourly  basis  during  the  contractual  period 
were  Miss  B.  Peeler  and  Miss  M.  Zenderman. 
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During  the  contractual  period,  four  Scientific  Reports  were 

issued: 

1,  Static  stability  measures  in  the  atmosphere,  by  W.  Lawrence  Gates. 
Scientific  Report  HO.  3,  31  October  i960,  22  pp.  (This  report 

was  prepared  under  the  previous  Contract  AP  19(6o4)-4965). 

Published  under  the  same  title  in  J.  Meteor. ,  18:  526-533 

(1961). 

2,  Preliminary  study  of  numerical  errors  in  the  integration  of 
barotropic  flow  on  a  spherical  grid,  by  W.  Lawrence  Gates  and 
Christopher  A.  Riegel.  Scientific  Report  Ho.  4,  15  March  1961, 

46  pp.  Published  under  the  same  title  in  J.  Geophys.  Res. 

67:  773-784  (1962). 

3,  The  propagation  of  error  in  numerical  Integration  of  the 
discrete  vorticity  equation,  by  W.  Lawrence  Gates.  Scientific 
Report  Ho.  6,  15  May  1962,  37  PP- 

The  effect  of  finite  differences  on  the  growth  rates  of 
unstable  waves  in  a  simple  barocllnic  model,  by  Christopher  A.  Riegel. 
Scientific  Report  Ho.  7,  1  October  1962,  24  pp. 


W.  Lawrence  Gates 
Project  Director 
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ABSTRACT 


Using  a  spherical  finite-difference  grid  system  covering 
the  hemisphere,  numerical  solutions  of  the  non-divergent  barotroplc 
vorticity  equation  have  been  obtained  for  periods  up  to  10  days  by 
relaxation  methods.  In  the  case  of  analytic  initial  conditions, 
the  solutions  are  sensitive  to  abrupt  changes  of  the  longitudinal 
mesh  size  in  the  regions  of  appreciable  amplitude  of  the  stream- 
function  tendency.  A  doubling  of  the  longitudinal  mesh  size  from 
5  to  10  deg  at  U5°N  produces  a  progressive  tearing  or  shearing  in 
the  stream function  which  is  quite  noticeable  at  10  days .  When  the 
grid  variation  is  removed  to  70°R,  however,  the  solutions  at 
10  days  are  free  of  this  defect,  although  exhibiting  the  expected 
phase  lag  relative  to  the  analytic  solution.  When  applied  to 
hemispherical  barotroplc  prediction  with  actual  500  mb  data,  this 
spherical  grid  scheme  appears  to  perform  as  satisfactorily  as  do 
conventional  (rectangular)  grid  formulations . 
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1.  Introduction. 

In  an  earlier  study  (Gates  and  Riegel,  1962)  the  problem 
of  numerically  integrating  the  barotropic  vorticity  equation  on 
a  spherical  finite -difference  grid  was  considered.  By  comparing 
the  numerical  solutions  with  analytic  solutions  for  a  simple  harmonic 
initial  condition,  the  dominant  sole  of  the  space  truncation  error 
was  identified.  It  is  the  purpose  of  the  present  paper  to  extend 
these  analytic -numerical  comparisons  to  somewhat  longer  periods 
of  integration,  and  to  report  on  an  improvement  in  the  grid  design 
as  a  result  of  these  investigations. 

It  is  sufficient  here  to  recall  that  the  barotropic  vorti¬ 
city  equation  for  horizontal  non-divergent  flow  may  be  written  in 
the  form 

2^5  =  _J _ r^s  ^  21  1 

dt  2<f>  ?4>  'JA  J  >  (1) 


with  the  vorticity  given  by 


_J _ 


5* 


) 


_J _ 


(2) 


Here  is  a  streamfunctlon,  r  the  radius  of  the  assumed  spherical 
earth,  ^  the  latitude,  and  A  the  longitude.  This  equation  is 
approximated  by  the  finite-difference  equation 


where 


ht  *  he  A*  ( )"‘  cos  <f>  (4) 

=  rn.  (2.)'1  sih^  ,  (5) 

and  where  A 4  and  are  the  basic  (5  deg)  latitudinal  and  longi¬ 
tudinal  mesh  increments,  respectively,  with  corresponding  integer 
indie  ies  i  and  j  . 

The  details  of  this  difference  scheme,  the  treatment  of  the 
pole,  and  the  method  of  numerical  solution  have  been  described 
earlier  (Gates  and  Riegel,  1962).  The  distinguishing  notion  of  this 
approach,  however,  is  the  identification  of  the  points  L±nt  in 
(3)  above  as  those  at  the  neighboring  points  ±  5  deg  long  away 
(i.e.,  m.  =  1)  for  0  s  ^  s  45N,  and  taking  m,=  2  for  the  latitudes 
50-65N,  hi  =  3  for  70N,  1H  =  4  for  75N,  m.  =  6  for  80N  and  m=  12 
at  85  N.  Equation  (3)  is  itself,  however,  applied  to  each  point 
of  the  basic  72  X  18  point  grid,  with  the  pole  treated  separately. 

In  this  fashion  comparable  grid  resolution  is  maintained  over  a  wide 
range  of  latitude  without  interpolation. 

2.  Extended  numerical  Integrations. 

In  the  previous  study  (Gates  and  Riegel,  1962)  it  was  found 
that  in  the  solution  of  (3)  by  the  extrapolated  Lietxnann  relaxation 
method,  an  overrelaxation  coefficient  4  =0.36  and  a  streamfunction 
tendency  residue  tolerance  of  0.07  percent  of  the  maximum  initial 
tendency  were  suitable  values  for  the  analytic  case  selected. 
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These  numerical  solutions  have  now  been  extended  to  10  days  (vith 
time  steps  At  =  1  hr),  and  are  shown  in  a  sixth  part  in  Figs.  1-5* 

A  sixth  part  of  the  associated  streamfunction  error  fields  are  shown 
in  Figs.  6-10  for  this  case  of  six  symmetrical  hemispheric  waves 
moving  eastward  at  20  deg  long  per  day  without  change  of  shape. 

It  may  be  noted  that  the  waves  retain  their  general  initial 
configuration  until  approximately  6  days,  and  move  at  about  96 
percent  of  the  true  wave  speed.  After  this  time,  however,  an  in¬ 
creasingly  serious  distortion  occurs,  such  that  by  10  days  the  initial  wave 
pattern  is  noticeably  altered.  Closer  examination  of  Figs.  6-8 
shows  a  steadily  increasing  streamfunction  error  field  of  the  type 
to  be  expected  from  a  persistent  phase  speed  lag  of  the  predicted 
versus  the  actual  wave.  Figs.  9  and  10,  however,  show  a  somewhat 
more  rapid  error  growth,  with  a  pronounced  error  asymmetry  by  10  days. 

This  distortion  appears  first  near  the  latitude  at  which  the  effec¬ 
tive  longitudinal  mesh  size  is  first  doubled,  i.e.,  between  V?N  and 
50N.  Due  to  the  Increased  truncation  error  the  predicted  wave  speed 
to  the  north  of  this  juncture  is  less  than  it  is  to  the  south,  and 
this  evidently  introduces  a  fictitious  "backward"  tilt  to  the  wave. 

This  effect  is  here  less  noticeable  at  the  higher  latitudes  (at  which 
similar  longitudinal  mesh  size  changes  occur )because  of  the  somewhat 
smaller  wave  amplitude  in  these  regions. 

This  characteristic  wave  distortion  also  occurred  in  the 
extended  integrations  with  a  doubled  (10  deg)  basic  grid  size,  as 
shown  in  Figs.  11  and  12.  The  larger  truncation  error  in  this  case, 
however,  evidently  masks  the  distortion  so  that  it  is  less  clearly 
displayed  at  10  days  than  in  the  case  of  the  5  deg  grid  discussed  above. 


Figure  1.  A  sixth  part  of  the  computed  stream  function  solution 
at  t  =  2  days  with  the  original  5  deg  grid,  with  Oi  =  O.36, 

2  -2  1 

€  -  0.50  m  sec  ,  and  At  =1  hr.  The  dashed  lines  are  the 

7  2  -1 

corresponding  analytic  solution.  The  units  are  10  m  sec  ,  and 
the  arrows  "A"  and  "C"  indicate  the  positions  of  the  analytic  and 
computed  stream  function  troughs,  respectively.  All  stream  function 
values  are  negative. 
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Figure  11.  A  sixth  part  of  the  computed  stream  function  solution 

at  t  =  10  days  with  the  original  10  deg  grid,  with  Gi.  =-  0.3 2, 

2  “2  72  X 

€  =  0.5r'  m  sec  ,  and  At  =  1  hr.  The  units  are  10  m  sec  . 

See  fig.  1. 


Figure  12.  A  sixth  part  of  the  error  field  of  the  computed  stream 
function  solution  at  "t/  =  10  days  with  the  10  deg  grid.  The  units 
are  lO^m^sec-^.  See  fig.  1. 
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With  the  10  deg  grid  (see  Gates  and  Riegel,  1962)  the  wave  speed 
in  the  numerical  solutions  is  lowered  to  approximately  80  percent 
of  true  speed,  and  is  further  evidence  of  the  dominant  role  of  spatial 
truncation  error  in  these  calculations. 

3.  Extended  numerical  integrations  with  a  modified  grid. 

In  an  effort  to  reduce  this  wave  distortion  in  middle  latitudes, 
the  differencing  schema  was  changed  such  that  in  (3)  -  (5)  the 
parameter  V*=  1  for  latitudes  0  £  ^  £  JOS,  2  for  7 5H,  h\  =  3 
for  80N,  and  m,  =  6  for  95N.  The  polar  calculations  still  use  only 
6  of  the  points  at  85IT,  fbrming  a  polar  hexagonal  grid.  A  sixth 
section  of  this  modified  differencing  grid  is  shown  in  Fig.  13, 
although  the  full  array  of  points  for  each  3  deg  lat  and  3  deg  long 
continue  to  be  carried  in  order  to  avoid  interpolation  (see  Gates 
and  Riegel,  1962).  The  values  of  the  coefficients  in  (3)  may  readily 
be  determined  from  those  given  earlier  for  the  unmodified  grid. 

Since  the  closest  points  now  consulted  in  an  application  of  (3), 
or  in  the  determination  of  ^  ,  are  approximately  190  km  apart  at 
70H,  the  value  of  At  was  reduced  to  £  hr  in  order  to  avoid 
computational  instability  with  the  conventional  centered  difference 
scheme. 

With  the  same  analytic  6 -wave  initial  conditions  used  earlier, 
a  series  of  tests  of  the  relative  efficiency  of  various  overrelaxat ion 
coefficients  was  made  for  the  previously  determined  streamfUnctlon 

O  mO 

tendency  residue  tolerance  €  =  0.30  a  sec  .  From  these  tests 
at  "t  =  0,  sunmarized  in  Fig.  14  along  with  the  corresponding  data 
for  the  unmodified  grid,  the  value  =  0.28  was  selected  for  use  in 
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Figure  13.  A  sixth  part  of  the  modified  5  deg  finite-difference  grid 
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the  extended  integration  on  the  modified  grid.  It  may  be  noted  that 
the  seemingly  small  modification  of  the  computing  scheme  has  had  a 
substantial  effect  upon  the  relative  efficiency  of  the  relaxation 
solution.  The  number  of  scans  of  the  grid  required  for  convergence 
after  the  Initial  solution  shovn  in  Fig.  14  was  approximately  15 . 

With  the  modifications  described  above  the  streamfunctlon 
integrations  to  =  10  days  are  shovn  in  Figs.  15-19,  and  the 
associated  errors  are  shovn  in  Figs.  20-24;  these  figures  correspond 
to  Figs.  1-10  for  the  original  (unmodified)  grid.  From  these  figures 
it  is  at  once  evident  that  the  numerical  solutions  obtained  vith  the 
modified  grid  are  superior:  by  10  days  no  vave  distortion  is  evident, 
and  the  relatively  smooth  error  distribution  reflects  the  systematic 
underestimate  of  the  vave  speed  due  to  the  dominant  spatial  trunca¬ 
tion  error. 

4.  Integration  vith  observed  atmospheric  data. 

The  use  of  the  spherical  difference  scheme  vith  analytic  data 
should  be  supplemented  by  an  Integration  using  observed  data  if  the 
scheme's  application  to  practical  numerical  weather  prediction  is 
to  be  examined.  To  this  end,  the  500  mb  contour  data  for  15  GMT 
18  November  1955  vas  selected  from  the  U.S.  Weather  Bureau's  Historical 
Map  Series,  and  vas  used  in  the  preparation  of  a  2-day  numerical 
integration.  The  observed  500  mb  contour  maps  for  this  period  are 
shovn  in  figs.  25-27. 

Before  the  integration  could  proceed,  the  streamfunctlon  ^ 
stipulated  by  the  assumption  of  horizontal  non -divergent  flow  first 
had  to  be  determined.  The  pertinent  equation  is  in  this  case  the 
balance  equation 
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Figure  13.  A  sixth  part  of  the  computed  stream  function  solution  at 

t  =  2  days  with  the  modified  5  deg  grid  (see  fig.  lj),  with  oC  =  0.28, 

2  —  2 

6  =  0.^0  in  sec”  ,  and  At  =  #  hr.  The  dashed  lines  are  the  corresponding 

7  2  -1 

analytic  solution.  The  units  are  10  m  sec  ,  and  the  arrows  "A"  and  "C” 
indicate  the  positions  of  the  analytic  and  computed  stream  function 


troughs,  respectively.  All  stream  function  values  are  negative. 


Same  as  fig.  15,  but  for  t  =  8  days. 


Figure  20.  A  sixth  part  of  the  e 


function  solution  at  "t  ®  2  days  < 
6  2 

unit 6  are  10  m  sec  .  See  fig.  r 
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■e  27.  The  observed  500  mb  topography  on  1500  GMT,  20  November  1955. 
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written  here  for  convenience  in  rectangular  coordinates .  In  a  number 
o t  studies  of  the  use  of  this  equation  at  500  mb,  the  non-linear 
terms  have  been  found  to  contribute  relatively  little  to  the  stream- 
function  solution  with  observed  contour  heights.  Jtor  the  present 
purposes  it  therefore  seems  sufficient  to  use  the  linearized 
version  of  (6),  which  in  spherical  coordinates  becomes 


_J _  2_ 


+ 


_j _ 

\r%Cos*$ 


1  4-i-^ 

J  9»*  H 


—  — j -  JL  (  CSS  4  —  )  ■+  1  _  ^ 

t^cos  ^  '  hacesx4  9AX 

The  introduction  of  centered  finite  differences  on  the  modified 

5  deg  grid  of  §  3  leads  to  the  difference  equation  approximating 

(7)  in  the  form 

(  h,  -Ha-t  *l3  )  (  +  (  *!,+•»*  -  H j  )  'fy- 1  +  *  i  (*{4  Mij  4  ^i-^-  ) 


-  (£hl+*hr,)+l*i  * 

+  \'  -  (*».+*<  >  _  <«> 

where  H,  and  ft  ^  are  given  by  (fc)  and  (5),  and  share 


*S  *  K,  i*  M  (*?j )-  ,  (9) 

with  the  Coriolis  parameter,  ^  the  Bossby  parameter,  and 
the  remaining  symbols  as  previously  defined.  At  the  pole  the  linear 


portion  of  (6)  is  approximated  on  a  locally  hexagonal  grid  as 

«ltvvVV*rwi  -  fr  (V W VV**' ),(10) 

where  the  subecripta  1  to  6  denote  6  equally-spaced  points  at 
85°!?,  and  the  subscript  P  denotes  the  pole. 

Starting  from  the  known  %  ,  equations  (8)  and  (10)  were 
solved  for  the  streaafunction  ^  by  extrapolated  Liebaaim  relaxation, 
using  t  -  0  as  a  boundary  condition  at  the  equator.  A  series  of 
tests  for  various  values  of  the  overrelaxation  coefficient  and 
residue  tolerance  were  made,  and  are  shown  In  fig.  28.  The  solution 
with  6  =  0.1  x  10^  m2  sec"1  and  0(  =0.29  was  ultimately  selected, 
although  there  was  very  little  variation  in  the  solutions  with  « 
for  this  value  of  6  .  We  may  note  the  similarity  of  the  useful 
(X  range  here  with  that  in  fig.  14,  and  conclude  that  the  h3 
term  in  (8)  has  caused  no  difficulty  in  the  relaxation  solution. 

The  relaxation  solutions  of  (3)  itself,  once  ^  had  been  found,  re¬ 
quired  approximately  110  scans  per  (£  hr)  time  step,  with  the  values 

2  -2 

oC  m  0.30  and  £  »  0*05  ®  sac  ;  these  constants  were  selected  from 
the  data  of  fig.  14,  although  €  was  here  lowered  by  a  factor  of 
10  in  the  interests  of  accuracy. 

In  order  to  verify  the  numerical  integrations  for  V)/  f  (8) 
and  (10)  were  rewritten  for  ?  in  terms  of  V|/  ,  and  again  solved 
by  the  relaxation  method.  In  these  inversions  the  overrelaxation 
coefficient  and  the  residue  tolerance  were  given  the  values  e<  =  0.29 
and  6  3  0.5  ft,  respectively,  and  required  approximately  590  scans 
for  convergence.  The  initially  observed  zonally  averaged  contour 
height  at  the  equator  was  used  as  an  equatorial  boundary  condition. 
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The  scans  of  the  1297-point  modified  5  deg  grid  required 


for  convergence  of  the  relaxation  solution  for  Vp  (actual  atmospheric 
500  mb  data  of  1500  GMT,  l8  November  1955). 
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The  forecast  500  nfc  contour  maps  computed  in  this  manner  from  the 
1-  and  2«<iay  numerical  integrations  of  (3)  with  the  modified  grid 
are  shown  in  figs.  29  and  30.  These  forecasts  are  not  of  high 
meteorological  accuracy  in  that  the  behavior  of  a  number  of  major 
wave  systems  over  the  hemisphere  is  not  well  predicted.  An  inde¬ 
pendent  integration  (not  shown),  with  the  same  initial  and  boundary 
conditions  on  a  conventional  rectangular  grid  system  of  comparable 
middle-latitude  resolution,  was  physically  inadequate  in  the  same 
instances,  however,  so  that  these  failings  may  be  attributed  to  the 
non-divergent  barotropic  model  rather  than  to  the  grid  scheme  itself. 

5.  Conclusions. 

In  using  a  spherical  grid  scheme  with  changing  longitudinal 
mesh  size  (Gates  and  Riegel,  1962),  it  appears  important  to  avoid 
an  abrupt  change  of  mesh  size  in  those  latitudes  at  which  large 
gradients  characteristically  occur.  The  latitude  70°K  appears 
suitable  for  a  first  doubling  of  a  5°  longitude  mesh  size  for  typi¬ 
cal  atmospheric  flow  at  500  mb,  and  this  spherical  grid  appears  to 
yield  a  smooth  and  well-behaved  solution  over  the  hemisphere.  It 
is  felt  that  such  a  scheme  may  prove  of  considerable  use  in  global 
numerical  weather  prediction  in  which  artificial  lateral  boundaries 
or  overlapping  rectangular  grids  are  not  desirable. 
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